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ABSTRACT 


The dependence (or independence) of baseline accuracies, obtained from 
a typical mission of a spacebome ranging system, on several factors is investi- 
gated. The emphasis is placed on a priori station information, but factors such 
as the elevation cut-off angle, the geometry of the network, the mean orbital 
height, and to a limited extent geopotential modeling are also examined. 

The results are obtained through simulations, but effort has been made 
to give some theoretical justification whenever possible. Guidelines for freeing 
the results from these dependencies are suggested for most of the factors. 
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1. INTRODUCTION 


In the past decade or so laser ranging to artificial satellites proved to be 
one of the most precise and efficient tools for geodetic positioning. The benefits 
from its use were realized early enough to encourage further development of the 
hardware as well as extensive application to problems related to earth dynamics. 
One of the areas where this system will be of major importance is plate 
tectonics . Almost half a century after Alfred Wegener published his continental 
drift theory (Wegener, 1928], space geodesy provided scientists with a sound tool 
for measuring the relative motion of the continental plates. As our knowledge and 
understanding of geophysical phenomena related to plate tectonics improved, it 
became apparent that there exists a high correlation between the location of the 
plate boundaries and earthquake epicenters. Further, it has been the conviction of 
several scientists that geophysical activity in the region of a fault contains vital 
information about the actual occurrence of earthquakes. It is therefore highly 
desirable to be able to monitor such activities (dilatancy, strain accumulation, tilt, 
etc.) as they are related to seismic hazards. The regional aspect of plate 
tectonics, therefore, is mainly concerned with the deformation of the plates near 
their boundaries at the fault zones. The best way of determining this deformation 
is monitoring the motion of several benchmarks located near the fault relative to 
that of points significantly away from it. Since a fault zone can be of quite large 
extent, in order to be able to deduce meaningful results a large number of points 
is required and quick, frequent rcsurveying of the area. 

A system that can meet all the requirements and still be cost effective is 
a Satellite Ranging System (SRS). The idea behind this system is the Inversion of 
the traditional satellite trilate ration scheme. Due to the large number of points whose 
positions must be determined, the active (and expensive) station, e.g., the laser, 
is placed on the spacecraft and the ground stations are targeted with relatively 
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cheap reflectors. The advantages of this scheme in terms of cost are obvious. 

Historical Review 

The idea of the "upside-down" laser was first suggested in the late '60's 
at the time NASA undertook the task of improving the existing hardware so that 
higher accuracies could be achieved. In 1974 a research team from The Ohio 
State University Department of Geodetic Science undertook jointly with the 
Smithsonian Astrophysical Observatory the investigation of what was then called the 
"Close-Grid Geodynamics Measurement System" (CLOGEOS). Thepart investigated 
at SAO pertained to systems and objectives, while the one at OSU to optimum system 
use. The points studied under the second part included station configurations, orbital 
configurations, observational accuracy and data reduction techniques [Mueller et al., 
1975], Following this, a much more realistic simulation study was published 
[Kumar, 1976] in which the variation of several of the aforementioned factors, as 
well as new ones, was examined in detail. It must be mentioned here that at that 
time there was no final decision taken either for the type of ranging instrument to 
be used or for the orbit of the carrying spacecraft. The variations in orbital 
configuration were based therefore on theoretical arguments and the selection made 
was otherwise arbitrary. As for the ranging system, it was implicitly assumed to 
be a pulsed laser without ruling out any other suitable candidates (e. g. , radio Doppler). 
Although these studies did not produce the final answers to the problems involved 
with the system, they clarified to a high degree most of them and set up guidelines 
for future investigations . 

On the other hand, SAO produced a final report [SAO, 1977] on the investiga- 
tions conducted by them pertaining to systems and objectives for CLOGEOS. The major 
findings of that study were the following: (1) Most promising systems for relative 
positioning at the 1 cm level are the pulsed laser and the radio Doppler, (2) System 
accuracy is hindered by insufficient knowledge of atmospheric effects and loss of 
information due to variable weather conditions. As a result of these two findings, 
further investigations on the above subjects were proposed. As for the objectives 
for the employment of such a system, it was emphasized that the main application 
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of the system should be the densification of geodetic networks in the fault regions 
and their subsequent resurveying in order to produce the required information 
leading to a four-dimensional (space and time) deformation model . Secondarily, 
the system could also be used for studying other phenomena of interest to geodesists, 
geophysicists, glaciologists, and related disciplines. 

The momentum acquired from these two investigations and those conducted 
independently at Goddard Space Flight Center [Vonbun et al. , 1975; Agreen and 
Smith, 1973] along with the recent developments in hardware capabilities pushed 
the investigation into the next phase. The use of a pulsed laser was finalized and 
the idea of testing the prototype using the space shuttle under development gave 
birth to a new system, the "Spacelab Geodynamics Ranging System" (SGRS). The 
San Andreas fault-system zone was selected as s test area and an error analysis of 
this specific system was performed at Business and Technological Systems, Inc. 
(BTS) under the guidelines set by GSFC [Gibbs and Haley, 1978], 

At this point it was felt that a variance analysis for the new system SGRS 
seemed proper. Our study was conducted in two phases. The initial phase is an 
analysis of the proposed experimental system with shuttle flights. The key issue 
in this analysis is the variation and dependence of the recovered baseline precision 
due to the variation of certain factors such as baseline length, a priori station 
information, observational accuracy, elevation eut-e'i' angle and network design. 

As this phase reached its end, a workshop at the University of Texas at Austin, 
organized by GSFC, brought together the various investigators of SGRS and the 
candidate users of the system. The purpose of this workshop was to review the 
current system design and gather information from its potential users pertinent to 
system improvement and operational system design. The results of the discussions 
and the recommendations from this meeting [Report from the Workshop on the 
Spacebornc Geodynamics Ranging System, 1979] proved to be of major importance 
for the design of the of the opt; rational system and for this reason a brief summary 
is given below: 

1. The i;.ser transmitter -receiver system must be designed to reduce the 
cost of the ground reflectors (below $1000 per unit) and still be capable of centimeter 
level geodesy. 
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2. Granted that the system is designed along the lines set during the work- 
shop sessions, geophysicists and seismologists concluded that several unique appli- 
cations of the system are possible and of great interest to the scientific world. 
Primarily, the system should be deployed at various areas around the globe in an 
attempt to "capture" a moderate-sized earthquake. 

3. Tectonic plate motion monitoring can be achieved by use of a system such 
as SRS both for near boundary deformations as well as for relative plate motion 
determination. 

4. In the light of the proposed system design, several— mainly of geophysical 
interest— experiments are proposed for the study of interplate and intraplate movements 
and their relationship to area seismicity. 

5. The system can be used to establish global and local geodetic control net- 
works of high quality. Mapping and resurveving of large areas with sparse or no 
geodetic control could be covered very rapidly, effectively and, above all, at low 
cost compared to classical methods. 

6. Various other applications of the system are possible (glaciology, atmos- 
pheric sciences, precise time transfer, etc. ) provided that the system be designed 
in a cost effective manner and be capable of achieving relative positioning accuracy 
of one to two centimeters over intersite distances of ten to fifty kilometers . Since 
the technology is available, it is recommended that the system be designed for a 
high flying (~ 1000 km mean altitude) dedicated satellite and equipped with a short 
pulse (0. 2 ns) laser. 

The analysis of the capabilities of the new system as it evolved from the 
guidelines set at the above meeting constitutes the second phase of the present study. 
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2 . OUTLINE O F THE INVESTIGATION 


The main objective of thir study is to determine the dependence of the 
recovered baseline precision on the following factors: 

(1) A priori station information, 

(2) observational accuracy, 

(3) geopotential model, 

(4) elevation cut-off angle, 

(5) baseline orientation (network geometry). 

The reason for selecting the baselines as representative end products of 
the whole process is twofold: primarily, the baseline lengths and the angles between 
them are the only estimable quantities in the adjustment and secondarily because in 
most applications of this system the conclusions will be based on the baseline length 
variation between missions. Since only a covariance analysis is performed, several 
simplifications were done in the course of simulating the observations. It must be 
pointed out that it was mainly due to restrictions imposed by the available software 
that we actually had to simulate observations (GKODYN requires either real or 
simulated observations in order to for m the normal equations and thereby compute 
variance-covariance matrices for the parameters or functions of the parameters 
such as the baselines). For a pure covariance analysis, no observations (real or 
simulated) are required. It was already mentioned in the Introduction that the 
pre sent investigation deals with two different versions of the basic system. The 
differences come mainly from the carrying vehicle and the selected orbit. In the 
first version it is assumed that the laser station will be placed in a low orbit 
(mean altitude ~ 100 km) aboard the shuttle (validation experiment). The second 
version is based on a free- flying dedicated satellite at a mean altitude of ~ 1000 km. 
The investigation of only these particular versions is based on the conclusions and 
recommendations of the SG RS Workshop f 19701 at Austin. 
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The characteristics of the active part of the system, the laser, were 
assumed to be the same in both cases. Detailed descriptions of the various 
components of the system appear in several reports compiled by the different 
agencies and companies involved in the development of the system. The only 
information which was actually used in this study is the rate at which the laser 
can operate, 10 pps, and the presently feasible resululier. for the single observation, 
10 cm. Although it is highly probable that the new generation pulsed lasers with 
0.2 ns pulsewidth and resolution of 2 cm will be operational at the time the 
system flies, we felt that it was more proper to perform our investigation 
based on present capabilities. The results of this study can be rather easily 
projected to indicate the impact of such an important improvement in the hardware. 

The simulation site selected covers part of the San Andreas fault system in 
California and Nevada. The 42 ground-based reflectors form a rectangular grid 
400 km long and 200 km wide (Fig. 1). The coordinates of these 'tations were 
obtained from BTS so that a direct comparison of our results with theirs would be 
feasible (Table 1). The observing sequence was the same for all passes and for both 
orbits. The sequence that we selected is shown in Mg. 2. In reality after a short 
acquisition period of about 10 s , the laser points to visible stations consecutively 
until all of them have been observed and then cycles back to make another set of 
observations. As the spacecraft ascends (or descends) over the horizon, only a few 
of the stations are visible. This results in an increased number of observations 
for the peripheral stations compared to the whole. It was felt, however, that for 
the purpose of this simulation and being consistent with other simplifications (e.g., 
weather effect simulation), the fixed schedule was adequate. 

The satellite orbits generated were based on a simplified force model. Two- 
body motion was assumed and |he gravitational earth model consisted simply of CM 
and ,12 • Onli ■ o secular variations due to .J 2 were considered in the numerically 
integrated equations of motion. The integration steps ize was 10 s. The simulated 
ranges were computed from interpolated state vectors using a third-order finite dif- 
ference method . This procedure was dic tated by the fact that in order to obtain state 
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Table 1 Participating Station Coordinate List 


Reference Ellipsoid: GRS 67 Earth-fixed Geocentric Coordinates 
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vectors at the 10 pps rate from the numerical integration, computational 
instabilities and unfavorable round-off accumulation introduced unacceptable 
biases on the results. The interpolation technique was tested against straight- 
forward integration and was found far more efficient in terms of computation 
time and accuracy. A summary of the information pertinent to the generation of 
the two orbits and the corresponding range observations is presented in Tables 2, 

3 and 4 . 

In order to be able to compare our results with those of other similar inves- 
tigations (like BTS), we used only nine passes (50% of the available) for the low orbit. 
The reason for using only one-half of the passes is that we assumed that there is a 
fifty-fifty chance that a pass will be observed depending on the weather conditions in 
the area. This, of course, is the simplest way of modeling the weather, but it was 
felt to be adequate for our purpose. For the high orbit only eight passes were used 
(25% of the available) and the observing rate was lowered to 1 pps . To compensate 
for this decrease in the amount of data, we increased the observational precision to 
10 cm//To sr 3.2 cm, assuming that the ten observations within each one-second 
interval are independent. This assumption had also been made at BTS when devel- 
oping their data set for the low orbit. We used this data set as provided by them for 
checking our software (GEODYN) and verifying their results. For the sake of brevity, 
the above three data sets will hereafter be referred to as OSUL, OSUH, and BTSL 
respectively. 

As mentio md previously, the computer pi'ogram used for the adjustment of 
the observations was GEODYN, obtained from NASA at an earlier time. This 
program is designed primarily for complex dynamic solutions and its use for the 
present investigation somewhat beats its purpose. However, since other investigators 
used this program already and a number of its capabilities simplified to a high degree 
the task undertaken, we decided to use it. We should mention here that in the present 
investigation only dynamic solutions were considered. Although in [Kumar, 1976] 
and [Kumar and Mueller, 1978] the geometric solutions are shown to be more promis- 
ing than the dynamic, we decided that on the basis of today's technology the realization 
of simultaneous ranging to the required more than six (or even four) ground stations 
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Table 2 Data used in the simulation procedure. 



Table 3 Distribution of Observations per Station per Pass for the 
Low Orbit (OSUL) 
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Table 4 Distribution of Observations per Station per Pass for 
the High Orbit (OSUH) 
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from a spacebome laser is practically impossible. Furthermore, we chose to use 
the short-arc mode, primarily in order to avoid the accumulation of biases in the 
recovered baselines due to uncertainties in the force field description (mainly the 
part pertinent to earth gravity modeling), and secondarily realizing that only a 
small fraction of a full revolution will be covered by observations since the stations 
are all spread over a very limited area. Unlike other software packages which are 
specifically designed for short-arc adjustment, GEODYN treats the short-arc solution 
as a collection of simultaneously reduced short (in duration) "long arcs." There are 
no approximate solutions to the equations of motion nor any other approximation 
whatsoever. As it is explained in [Mueller et al. , 1975) , one cannot expect to 
determine anything but relative positions from such a limited mission. The geo- 
potential model therefore must be held fixed and the optimum way for its description 
must be determined through suitable experiments. For our purposes we nominally 
used the GEM7 spherical harmonic expansion up to degree and order sixteen (16, 16). 
This was dictated by the fact that the same model was also used by BTS in their 
investigation. The groundtracks of the generated satellite passes are depicted in 
Fig. 3 (OSUL) and Fig. 4 (OSUH) . 
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3. THE ESTIMATION PROCESS 

The estimation process is described in [Martin et al . , 1976], This discus- 
sion, though, is limited in theoretical details and treats mainly the technicalities 
and the implementation of the process in the associated computer program. A simi- 
lar brief recollection of the formulas is also given in [Gibbs and Haley, 1978], Ac- 
cording to these references the method employed in GEODYN is Bayesian estimation. 
This type of estimation makes use of a priori information about the parameters in the 
form of a weight matrix derived from the prior distributions of the parameters. 

In most cases no such distributions are available and the weight matrices are 
derived from the assumed standard deviations of the a priori estimates for the 
parameters. In very rare circumstances a full variance-covariance matrix 
obtained from a previous solution is used for the determination of these weights. 

This latter procedure led several scientists into a different interpretation of the 
process, namely the so-called "least squares estimation with 'observed' parameters. " 
By doing this, assuming that our approximate parameters are the outcome of 
some measurements, we effectively change their character from fixed quantities to 
stochastic ones. It is not that there exist no problems where the parameters are 
inherently of that nature, but in our problems, especially when dealing with parameters 
such as Cartesian coordiaates, we cannot justify such an assumption. We should 
also point out that quite frequently the above procedure is used to improve 
the condition of the normal equations and in extreme cases to produce a full rank 
normal matrix for a problem which otherwise would be unsolvable in the domain 
of Cayleian matrix algeb>-j . The application of such "weighted" constraints on 
the estimates may very well distort the results in various ways and sometimes it 
may even result in unacceptable answers. From this point of view, the estimation 
process is directly related to the concept of estimable parameters, due to R.C. Bose. 
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It is proved in [Rao, 1973, p. 224] that the necessary and sufficient condition for a 
linear function of the parameters to be estimable is that the rank of the normal equa- 
tions is equal to the number of the parameters under estimation. A necessary and 
sufficient condition for the case where generalized inverses are used is also given in 
the above reference. The advantage of dealing with estimable quantities is that they 
are unique and unbiased for any solution of the normal equations and have minimum 
variance among all linear unbiased estimates. In this sense it is obvious that if 
we can identify the estimable parameters in the problem and if we can modify the 
model so that all the parameters are estimable, we have guaranteed ourselves a 
unique minimum variance unbiased estimate. This being the case, the need for 
a-priori statistical information on the parameters is alleviated, and our estimates 
are based only on the information provided through the observations. 

In the case where we cannot find such a set of estimable parameters to 
describe the system, we can still by-pass the need for exterior information by 
use of a generalized inverse solution. If it is possible to find linear parametric 
functions that fulfill the estimability criteria, then the analysis of the system can 
be done on the basis of these results. In the case of a geometric solution in 
satellite goedesv, for example, it has been shown [Mueller at ->1, , 1975] that the 
baseline lengths in the network and the angles formed by any Jvee stations are 
estimable. For most of the applications in geodynamics, these two quantities are suf- 
ficient for inferring motions and deformations in the area. The advantages of 
parametrizing a problem in terms of estimable quantities have been recognized 
by most scientists and effort has been made recently to identify these quantities for 
various geodetic problems. The complex functional relationships between the 
observables and the parameters are the only reason that estimable parametrization 
has not been in wide use yet. The case of laser observations, however, has been 
treated extensively in [Van Gelder, 1978], and the estimable parameters have been 
determined for extremely simplified models as well as for more complicated ones 
where secular perturbations due to J 2 were included. 

One therefore has all the tools to perform a proper study of a system such 
as the SRS, at least for simple simulations as is the case here. 
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If this is so, then onecan naturally question the reason behind this investigation 
and its major concern for the effects of a priori Information for nonestimable 
parameters. The answer to this is rather simple. Re?*’, world problems are 
far more complicated than an error analysis based or. a simulation of the system, 
and although major effort is currently devoted to improving and updating our 
estimation procedures we are far from being able to solve our problems in the 
fashion described above. Our best alternative, therefore, is to study in detail the 
current procedures to a degree that would allow us to justify our results and to 
determine how optimistic or pessimistic they may be due only to biases introduced 
by the employed procedure. 

In an effort to examine and clarify the inherent characteristics of Bayesian 
ef ’ imatlon, some of the major variations of this process are presented and com- 
pared to the well-known least squares estimation in the Appendix. It should be 
pointed out that the use of Bayesian estimation over least squares is an open 
question for statisticians as well as for scientists using these methods. Subjectivity 
versus objectivity is a rather philosophical question, and we feel that it is not the 
purpose of this study to provide the answer. One conclusion that can be drawn is 
that there are situations in applied science where one method is better than another. 
It is therefore our responsibility to choose between the two and to do this we must 
study both. The stand we take here is that the choice will be made on the basis of 
the problem requirements only, irrespective of the personal preference of the 
investigators. 


4. CONSTRAINTS, RANK DEFICIENCY AND ILL-CONDITIONING 
IN SHORT -ARC SOLUTIONS 



In the course of this investigation we have pointed out that we arc primarily 
Interested in determining the influence of station- related a priori information on 
the baselines' precision, it is only natural, however, to address the following 
question: What is the rationale for using such information? The answer to this 
question is not as simple as one might expect. If the problem is studied in depth, 
it will soon become apparent that this is just another way of posing the following 
fundamental question: Subjective (Bayes) or objective (Gauss- Markov) estimation? 
This is an open question for the statisticians and has a rather philosophical than 
mathematical nature. One should, therefore, not expect an answer from a limited 
study as this one. We would rather discuss the geodetic aspects of the problem and 
refrain from attempting to provide a conclusive answer, as that is outside the scope 
of this study . 

In the theory of linear spaces, rank defect of a matrix is the number of linear 
dependencies which exist among its columns. If we consider a matrix A as a linear 
transformation from a space JR* to a space ■ \ that is, A : M* R r , then the 
following example illustrates the concept of rank deficiency. From linear algebra 
it is known that the column rank of a matrix is the same as its row rank, where 
by column (row) rank we mean the maximum number of linearly independent columns 
( o -. >) of the matrix. Suppose n m for A and rank (A) k where the dimensions of 
A arc n by m (rows x columns). From the above theorem it is obvious that the rank 
of A will be at most m, when all its columns are linearly independent. If, therefore, 
k - m, the rank deficiency of A is zero or A is of "full rank." If k < m, then the 
difference m - k m - rank (A) is the rank deficiency of A. In terms of linear 
operator theory this means that the null space of the operator represented by A is 
not just the zero element, Ker(A) / 1 0 j, where the null space or kernel of a linear 
operator on ■* is defined as the subset of its domain M*, consisting of elements x / 0 
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for which Ax - 0 holds. Using another theorem new, we can stale that the trans- 
formation is not "one-to-one" or "tniectlve" and the equation Ax y therefore does 
not have a unique solution for x. 


The last remark provides the link between the mathematical and the statistical 
and physical interpretation of rank deficiency. If y denotes the vector of observations 
in a system and x the vector of para neters, with A the design matrix relating the two, 
then given a set of n independent observations (n ' m), only k rank(A) parameters 
out of the total m can be estimited uniquely. Physically this simply means that the 
given set of observations does not contain the information needed for the determination 
of the m - k parameters. In this sense, the parameters which can be estimated 
from the given observations constitute the set of "estimable" parameters of the 
system. We can conclude therefore that for each kind of observation in a given 
system there is a corresponding set of "estimable" parameters. It is easy to see, 
for example, that observing the velocity of a vehicle is not enough to define its 
position in space and time. We must find, therefore, a way to remedy this handicap 
in order to solve the problem. This isbasically done through the application oi 
constraints on the nonostimablc parameters. We can eliminate, for instance, these 
parameters In adopting certain values for them (absolute constraints), determined 
from a different approach. In this sense the "none stimn bio" parameters are not 
parameters anymore but constants of the problem. A solution of this problem 
through least squares (to account for redundant observations when n ■ m) will 
provide minimum variance unbiased estimates for ail (estimable by now) parameteis 
[Jiao, ld7d| . 

There are situations, however, where there is reason to believe that the 
available statistical information alxnit the nonostimablc parameters is of poor quality 
and enforcing them in the estimat ion process would result in unreasonable distortion 
of the results and sometimes even indicate nonexistent inconsistency In the obser- 
vations. The solution of the problem in this ease can lx 1 obtained utilizing what 
is called a set of flexible constraints. These in turn can be either arbitrary in 
numlier (but enough to provide a solution) or they can lx- just enough to make the 
least squares normal equation matrix (N A A) of full rank. The* latter are 
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called "minimum" or "minimal" constraints. It is obvious that these minimum 
constraints are not unique and the stability of the solution depends on the selection 
of these constraints. There is, however, a set of minimum constraints which 
is best in the sense that it does not discriminate among the parameters and provides 
a solution which is least affected by numerical instabilities in the system. This set 
of constraints is called "inner adjustment" constraints, and the theoretical and 
practical aspects of their application have been discussed extensively in [Rao,1973] 
and [Blaha, 1971]. We will not attempt a discussion of the pertinent details again, 
but we will point out that this approach is basically the same as solving the original 
set of normal equations by use of a generalized inverse of N [Rao, 1973]. The 
advantage of the "inner constraints" is that we obtain the same results without the 
troublesome computation of a generalized inverse. The application of this type of 
constraint in the present investigation was not possible due to limitations in the 
available software. 

Finally, another type of constraint which is widely used for the solution of 
rank deficient systems is the "weighted" or "relative" constraint. Since these 
were used in our test, we will elaborate on them and give some more details in 
addition to those which were already mentioned in the discussion of the estimation 
process. The idea of "weighted" constraints originated from the Bayesian approach 
of estimation in linear models . In this case, however, the a priori information 
which is added to the normal equations is based on an a priori known distribution 
of the parameters and the reason for using this information is well justified if 
one accepts Bayes' theorem. What we are trying to stress here is that in Bayesian 
estimation, the weighting of the parameters on the basis of a priori information is 
not intended to alleviate the rank deficiency in the problem nor the ill-ccnditioning 
of the normal equations. It merely makes use of all available information in order 
to arrive at estimates which are closer to the true values of the parameters, 
although not unbiased anymore. If, therefore, we know the prior distribution of 
the parameters, then this approach is fully justified as long as we are aware of the 
consequences of the Bayesian approach (biased parameters) . In this sense the 
inclusion of a priori information on the parameters should not at all be related to 
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the inherent rank deficiency or ill-condition of the problem. 

Weighted constraints, however, are used to overcome this problem. Their 
application in this sense is very dangerous as far as their effects on the results are 
concerned. In general, the argument on which this practice is based is that prior 
knowledge of the parameters' range of values is available through direct observations 
on them or from a previous solution. In the first case we change the role of the 
parameters to observations, and by doing this we are effectively removing them from 
the parameter list altogether. If in such a case we applied the weighted constraints on 
the nonestimable parameters only, then we have resolved the rank deficiency of the 
problem since we already know that for estimable parameters the design matrix 
(hence the normal equations) is in general of full rank. The catch is that in several 
cases we are either unable to observe the parameters because of their nature (e.g. , 
Cartesian coordinates) or we directly substitute their "observed" values with some 
approximate ones. In both these cases, the variance of these "measurements" is 
based on some personal confidence interval for the assigned values rather than 
what actual measurements would indicate. In the limit, as the a priori variance 
is decreased, the results of this adjustment are equivalent to those of absolute 
constraints, where we have changed the role of the nonestimable parameters to 
some adopted constants of the problem. For the case where we indeed have direct 
observations on some parameters (e.g., absolute gravity measurements in a gravity 
network), then we must treat them as observations with the proper variance- 
covariance matrix as we do for all other observations . This, of course, brings in 
the problem of determining the relative weights when different types of observations 
are simultaneously adjusted. This is a different problem, however, and we need 
not concern ourselves with it at this point. The important thing to note is that 
this type of observation should be treated in the usual manner and not used as a 
means for circumventing the rank deficiency in the problem. 

As far as the second case is concerned, where the weighting is based on a 
previous solution, we can identify two subcases. If all or some of the parameters 
were obtained from a previous adjustment, then their full variance-covariance 
matrix should be used in the new one, and the result is a standard Bayesian adjustment 
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as described previously. If, however, only the diagonal elements (variances) are 
used (which is common practice), then the resulting parameters (which are biased 
due to the use of a priori information) are not minimum variance - minimum bias 
estimates since the arbitrary diagonalization of the covariance matrix has changed 
their prior distribution from the one represented by the full matrix. Such practice 
is dangerous not only for the rank deficient problems but for those with full rank 
as well. In the case of the first, it merely provides deceiving results while in 
the second case it may distort the results instead of improve them. Some investi- 
gators choose to ignore the consequence of the last remark since they only judge the 
quality of their results from the magnitude of their a posteriori variances. These 
variances, however, are affected (and usually deflated) by the wrong a priori inputs 
and their statistical interpretation is very much questionable. 

This discussion can be summarized in the following. Rank deficiency is 
an inherent characteristic of each observational system, and in order to over- 
come it we must study the system and determine the source of this deficiency. 

This can be done through the determination of the estimable parameters for the 
system and their comparison with the list of our soive-for parameters. If there 
are observations (direct or indirect) which can be done to determine the non- 
estimable parameters, then we can perform these observations and include them 
in our adjustment alleviating the rank deficiency. When this is not possible but 
we have prior information in the form of a full covariance matrix, then we can 
perform a minimum variance - minimum bias estimation (Bayesian approach) 
where we obtain a solution for all parameters at the expense of their unbiased- 
ness. Alternatives to these are the inner constraint approach or a direct gen- 
eralized inverse solution and the adoption of absolute constraints. We point out, 
however, that the first does not provide unbiased estimates for nonestimable 
parameters. The similarities of this approach to Bayesian estimation are 
pointed out in the Appendix. An interesting and illustrative discussion of this 
approach is given in [Grafarend and Schaffrin, 1974]. The implications of the 
use of absolute constraints are rather obvious. The results are as good as the 
adopted values. One only needs to consider the by now popular leveling network 
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example where the absolute heights of the points can take different values 
depending on a single point's adopted height [Van Gelder, 1978]. 

Having discussed the nature of the various approaches to overcome rank 
deficiency problems in general, we proceed in examining the work and the results 
obtained therefrom for a satellite laser ranging system, in particular when used 
in the short-arc mode. The theoretical investigation of this system is given in 
[Van Gelder, 1978] where it is shown that for the simple case of a general elliptic 
orbit secularly perturbed by the J 2 harmonic, the rank deficiency is two. From 
the numerical tests which were performed for the case of the spaceborne laser, 
however, the results seem to contradict theory. Test runs with even three absolute 
constraints on a single station's position vector failed to provide an acceptable solution 
indicating serious instabilities in the normal equations. This, however, should not be 
taken as a proof that there is a flaw in the theoretical investigations, but rather as an 
indication that for a strong solution there is something more than rank deficiency to 
be considered. That other element was already mentioned in the estimation process 
discussion, and it is related to the geometry of the problem. In our investigations 
we are dealing with an extremely limited area (only 200 km by 400 km) and with 
short arcs of length which at best reach only a tenth of a full revolution. With such 
poor geometry any information about the coordinate system definition, which would 
normally come from the orbital dynamics, is so little and insufficient that the 
ill- conditioning of the normal equations dominates the problem rather than the 
inherent rank deficiency. In this case it may be that by increasing the number of 
observed arcs this instability is greatly reduced. This, however, must be investi- 
gated since such an increase would also introduce new parameters in the original 
problem. At present a definition of the origin and orientation of the coordinate 
system through six suitable absolute constraints seems to be our best alternative. 

Since these are a set of "over-constraints" (the theoretical rank deficiency still 
remains two) but essential in order to provide results clearly independent of 
numerical instabilities, we refer to them as "quasi-minimum" constraints 
following [Van Gelder, 1978], 
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5. DESCRIPTION OF THE EXPERIMENTS 


The investigation was conducted in three stages which offers a natural 
classification for the description of the experiments. In the initial stage preliminary 
experiments were conducted in order to familiarize ourselves with the problem, the 
available software and to set up guidelines for experiments that would follow. The 
second stage deals with the experiments on the low orbit and the final stage with 
similar experiments for die high orbit. 

A, Preliminary Experiments 

When a simulation study is conducted, it is very important to have the ability 
to check the results with those obtained independently either through a different 
approach or a different study by another organization. Since in our case a similar 
study was conducted in parallel at GSFC and BTS, we followed their guidelines in 
setting up the experiments in order to be able to compare our results. The main 
purpose of the preliminary experiments was to become familiar with the problem 
and to test our own software (GEODYN version 7508.0) using the data provided by 
BTS (BTSL). Since our main concern is the effect of a priori information on the 
baseline precision, the preliminary test focused on the variation of this factor. 

These tests explored the baseline precision variations as the a priori information 
about the stations, the orbit and the observational accuracy varied (Table 5). 

The effect of an erroneous geopotential model on the baseline precision was also 
investigated. As is pointed out in [Gibbs and Haley, 1978], the BTSL data set 
was generated on the basis of GEM1 (4,1), while in our tests the earth is 
modeled through GEM7 (16, 16). We felt that this inconsistency between the data 
generation and the reduction techniques should be cleared as to its effect 
on the precision of the baselines. Intuitively one expects that since we are dealing 
with a limited area and we are doing a short-arc adjustment, there should be no 
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Table 5 Baseline Standard Deviations! Preliminary Test Results. 


Test 

! OSU 1 

OSU 2 

OSU 3 

OSU 4 

OSU 5 

BTS la 

Station Constraints 2 

1 m 

1 m 

5 m 

25 m 

Quasi-Minimum 

Constraints 

1 m 

Orbital Constraints 3 

Constrained 

Constrained 

Constrained 

Free 

Free 

Constra ined 

Standard Deviation 
^*^»4jfObsorvatlons 

10 cm 

100 cm 





Baseline/ Length 

10 cm 

10 cm 

10 cm 

10 cm 

91-81 (26 km) 

0.98 

9.42 

0.98 

0,98 

0.98 

0.97 

91-71 (il km) 

1.01 

9.71 

1.01 

1.01 

0.95 

1.00 

91-61 (100 km) 

1 04 

10. 12 

1.04 

1.08 

0.98 

1.05 

91-51 <:!02 km) 

1. 14 

10.75 

1.14 

1.17 

1.04 

1.12 

91-21 (302 km) 

1. 17 

11. 10 

1.17 

1.30 

1.08 

1.16 

91-11 (403 km) 

1. 14 

10.06 

1.14 

1.33 

0.76 

1.11 

91-15 (449 km) 

1. 17 

10.09 

1.17 

1.30 

1.01 

1.17 

91-92 (52 km) 

1.08 

10.34 

1.08 

1.08 

1.11 

1.07 

91-93 (103 km) 

1. 04 

10. 15 

1.04 

1.08 

1.04 

1.06 

91-94 (152 km) 

1. 08 

10.28 

1.08 

1.11 

1 .08 

1.08 

91-95 (203 km) 

1.04 

10.24 

1.04 

1.11 

1.04 

1.06 

Average A Posteriori 







Standard Deviations of 

~ 20 cm 

~ 30 cm 

~ 80 cm 

~ 500 cm 

~ 4 cm 


Station Coordinates (X, Y,Z) 








l Table values are In centimeters. 2 Standard deviations in each Cartesian coordinate. 
3 Constraints as applied by BTS: 58 m In X,Y,Z and 5.8 cm/s in X, V, Z. 






change in our precision estimates. Indeed our numerical tests showed no 
difference (at least in the order of millimeters) in the baseline precisions 
although there are appreciable changes in their lengths. The results are not 
tabulated since they are identical for all baselines. 

We can summarize the findings of these preliminary tests in the following. 
An almost linear relationship exists between the precision of the observations and 
the precision of the recovered baselines. Although this is not a peculiarity of 
dynamic solutions (as it is for the geometrical ones) we can probably attribute it 
to the fact that there is a uniform distribution of stations in the area with approxi- 
mately equal numbers of observations from each station to an optimal set of 
satellite passes covering it as an umbrella from all possible angles. This whole 
setup strengthens the geometry of the problem to such a degree that in this respect 
it behaves almost as a geometric problem. As can be seen from Table 5, there is 
some variation in the baseline precision as the constraints were varied from case 
to case. In connection with the estimability problem in the short-arc mode, the 
last entry in the table shows an average standard deviation of the recovered 
Cartesian coordinates of the stitions. It is quite obvious that since they are 
nonestimable quantities their precision is strongly dictated by the input a priori 
information. These numbers were not included to prove that these quantities are 
nonestimable, but rather to show the pitfalls one faces if he were to choose 
these quantities as the basis of his investigation. The a priori weighted constraints 
on the orbital elements were introduced in accordance with the guidelines and 
tests conducted at BTS [Gibbs and Haley, 1978]. For reasons discussed in 
[Van Gelder, 1978], these constraints were never again introduced in the system 
in the tests that followed. 

In the last column we included the results obtained at BTS in their Run / la 
[Gibbs and Haley, 1978, p . Gf> | which most closely resembles our Test # 1. The 
tabulated values in the above reference pertain to precision of the horizontal com- 
ponents of the baselines only, but for comparison purposes and because of the 
relatively short distances involved they can be taken as the precisions of the base- 
lines' lengths themselves. It is evident from these numbers that our results are 
in excellent agreement with theirs. 
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B. The Low Orbit Experiments 


The results of the preliminary tests indicated that most of the variation in 
the baseline precisions came from the variation of the weighting scheme and the 
selection of an elevation cut-off angle. The dependence of the results on the 
observational precision turned out to be rather straightforward, and therefore 
no farther investigation for this factor seemed necessary. As for the network- 
orbit configuration effects, related to the geometric strength of the problem, 
little could be done since the orbit and the network design were given and assumed 
to be the common denominators for all experiments. 

With the above in mind the investigation concentrated mainly on the effects 
of weighting and elevation cut-off angle and to a lesser degree on the effects of 
geopotential model variations. Three different weighting schemes were adopted: 

(1) strong weighting for all stations, assuming that the standard deviation for each 
of the three coordinates {X, Y, Z) is 1 m, (2) mild weighting based on a standard 
deviation of 25 m in each component, and (3) "quasi-minimum constraints" solutions, 
that is, the coordinate system is defined by holding fixed six coordinates properly 
distributed among three selected stations . 

The "quasi -minimum" constraints for both OSUL and OSUH orbits were 
arranged as follows: 

(i) Station REF011: X and Y coordinates held fixed, 

(ii) Station REF015: Y and Z coordinates held fixed, and 

(Hi) Station RE F09 5: X and Z coordinates held fixed. 

The elevation cut-off angle was varied twice, 20° assumed to be the nominal value 
and 35° above horizon. The attention given to this factor is warranted by the fact 
that a higher cut-off angle could simplify the design of the ground reflectors and 
reduces their manufacturing cost. Some of the teres were repeated using a different 
geopotential model in order to verify the results obtained from the preliminary 
experiments and detect possible interactions due to the variation of the other two 
factors. In all cases no a priori information on the orbit was introduced. Addi- 
tional discussion for the "quasi-minimum" constraints is left for a later section. 
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C. The High Orbit Experiments 


The results obtained from the investigation of the low orbit indicated that 
the designed experiments were suitable enough to allow the determination of the 
effects of considered factors on the baseline precision. In addition to this, the 
fact that a comparison of the results from the two orbits would reveal their 
dependence on the spacecraft's mean altitude convinced us to follow the same 
experimental design. The description of these experiments is omitted since they 
are identical to those for the low orbit as described in the previous section. Addi- 
tional tests were performed in this case where the main objective was the determina- 
tion of the "effective" rank deficiency of the problem. By "effective" rank deficiency 
we mean the combined effect of the inherent rank deficiency of the short-arc mode 
adjustment and the deficiency arising from the ill-conditioned normal equations 
which is characteristic of the specific problem under study. The strategy followed 
in this case was to relax the number of ’'quas l -minimum" constraints and then 
examine the condition of the resulting normal equations. As one can gather from 
the above, the effective rank deficiency depends grossly on the design of the experi- 
ment and the quantity, structure, and quality of the collected observations. It would 
be useless, therefore, to try to quote for it a number such as two, three or six, 
while it is more appropriate to give some general guidelines that can be followed 
for a broad class of problems . 
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6 . RESULTS O F THE LOW AND HIGH ORBIT EXPERIMENTS 

The results which are presented and discussed in the following sections 
were obtained from the experiments conducted with the low orbit (OSUL) and the 
high orbit (OSUH) in the last two stages of the investigation. As explained in the 
previous section, there is a total of twelve basic experiments, six for each orbit 
(three weighting schemes for each of the two cut-off angles). The quantities ana- 
lyzed are naturally the standard deviations of the estimated baseline lengths. For 
a network of say n stations, there is a maximum of m - baselines which 

can be formed in all possible combinations (without repetition) among the stations. 
In our case n = 42, which yields m = 861; and considering the fact that there were 
twelve different solutions, we arrive at the total number of standard deviations to 
be analyzed: 10,332. Although theoretically a large data sample has several 
advantages, practically one can infer almost the same amount of information (not 
as accurately though) by restricting the analysis to a smaller sample formed on 
the basis of certain justifiable assumptions. We will elaborate on this a little 
further in the course of explaining the method of analysis, since these assumptions 
lay part of the foundation of our conclusions. 

A. Detection of Sources of Variation in the Sample Through 
an Analysis of Variance (ANOVA) 

The analysis of variance (hereafter referred to as ANOVA ) is a statistical 
method of analyzing measurements depending on various kinds of effects (called 
factors ) which simultaneously affect them, in order to make qualitative and quanti- 
tative Inferences of these effects (Scheffe, 1959). 

Since the method implies the existence of some measurements, it is only 
natural to expect that the first step— the setup of the experiment— poses an experi- 
mental design problem. The fact that inferences are to be made on the basis of the 
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experimental results leads in turn to a problem in estimation and decision theory. 

For a successful experiment one should identify the factors affecting the measure- 
ments and make sure that the designed experiment will yield all possible combina- 
tion treatments . Once the measurements are available they are arranged in a 
rectangular array which is basically a pictorial summary description of the experi- 
ment and the measurement process. Usually the factors in an experiment will be 
varied within a certain range of values which are of interest to the person who con- 
ducts the experiment. These variations constitute the levels for each factor. To 
obtain meaningful results we must have several observations (and even better, an 
equal number) at each level for all possible combinations. When more than one 
factors enter the problem, a setup such as that described above is called a 
complete factorial experiment . An example of such an experiment with three factors 
A, B, and C at I, J, K levels respectively is shown in Table 6. The simplest entity 
in an ANOVA table is the observation. Each observation is indexed in the following 
manner: one index for each factor plus one index which denotes the order of the 
observation within the M observations performed at each level combination (we 
tacitly assumed an equal number, M, of such observations in all levels). The next 
entity which is of interest to us is the set of the M observations for each treatment. 
They are easily identified in the table since they all have the same index values for 
the indices associated with the factors. This entity is called a cell . In the simplest 
setup there will be only one observation in each cell, (M ~1). 

There are two ways of Interpreting the ANOVA table. In the most usual case 
we assume that given R cells with M observations in each one, each cell represents a 
random sample of size M drawn from R normal populations with identical variance 
(7 2 . This means that we are only Interested in the specific variations (levels) of the 
factors as entered in the experiment disregarding the fact that these variations may 
be only a subset of many more possible. The effects of these factors on the observa- 
tions are therefore fixed with respect to the specific experimental setup, thereby the 
name jf this ANOVA model: fixed effects or Model I ANOVA. In the second case 
the interpretation is similar to the above except that now we make the assumption 
that the selected levels constitute a randomly selected sample from a large normal 
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Table 6 Arrangement of Data for a Three- Factor Experiment 

Factor A: I levels; Factor B: J levels; Factor C: K levels, 
M observations per cell (treatment or level combination). 



1 Levels of A 

L .. i 

• • • 

i 

Levels of B 

D 

II 

■a 

D 

l 

U 

j 

Levels 

of 

C 

1 

yim 

ym- 


yuu 
yu in 


ynu 

ym- 

i 

yuu 

ym- 

■ 


. 






K 

yiwi 

yn*M 

• • • 

yu«i 
yu* »* 


yn*i 

yu** 

i 

yum 

yu*- 


population of realizations of the considered factors, with an associated variance 
or? . As it is probably obvious by now, the names associated with this interpretation 
are: Random effects or Model II ANOVA. For the sake of completeness we mention 
here that it is possible to have an experiment where both fixed and random effects 
may be present simultaneously. This setup is usually referred to as the Mixed Model 
ANOVA. The next step, irrespective of which model we are using, is the estimation 
problem. Since this part is mainly computational we will restrict ourselves to 
discussing the procedure which should be followed for the example shown in Tabic fi. 
A generalization of the method for factorial designs of higher order than three will 
then be obvious without need for a general and involved presentation. 

The estimation process is based on the Gauss-Markoff theorem and it involves 
the computation of a number of "sums of squares" (SS) of deviations about various 
pivot values which will be explained in the following. These sums of squares will be 
used then in the computation of the "mean squares" (MS) which form the input of the 
decision theory problem. It can be easily verified that for the setup of Table 6 the 
total number of observations is given by the product IJKM, and the given set of data 
is said to span an UKM-dimensional space. Denote by { yi Jk » } the set of observa- 
tions, yun being the m^ 1 observation in the l,j,k "treatment combination" or the 
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cell with "coordinates” i,J,k in Table 6. Let Ut» be the expected value of the 
measurement on this combination of factors and iei* a ) ~ N (0, o 2 ) (independently) 
the random noise component of yijk«. Our linear hypothesis, which constitutes the 
basis of the estimation is: 

I yijk« * u t jk + eijkt 

n. ' 

I lej jk« 3 ~ independently N(0, a 2 ) 

Under ft each ui Jk is modeled as follows: 

ut j b = M 4 at + aj + a k + ajj 4 aj k + a ll( + atjk 

where M is the "grand mean" of the population, that is, the average of all the cells' 
u's . In ordc;- to avoid cumbersome notation involving multiple summations, aver- 
aging over a subscript is denoted by replacing the corresponding tndex by a dot. 

For example, = u and u tJk - y <Jk . . It turns out that under ft, an unbiased 
estimate of M is the average of all the observations, i.e. , M - y. . . . The 
remaining components in the expression for Usj k are the effects of the factors 
(a* , a*, a') and their "interactions" (ai 3 ®, a* k , ...» ai» c ), the effects, that is, 
which exist due to the simultaneous operation of these factors. The interactions 
i cannot be attributed to a single factor but only to their coexistence. If no inter- 

^ actions exist, then the factors are said to be "additive." In a sense there is no 

correlation among them and each operates independently of the others. When 
the number of observations in each cell is the same, as is the case here, the 
observational spaced can be decomposed into 2 r + 1 (provided M > 1) mutually 
orthogonal subspaces with a consequent decomposition of the total sum of squares 
into simpler sums . The dimensions of these subspaccs are the degrees of freedom 
(DF) associated with the corresponding SS, The estimates of the various effects 
and interactions are obtained from combinations of averages over different indices 
depending on which component we are estimating. For example. 


I 

l 

s 


A = y. ... 

at = yt... ~ y 

at 3* = yi 3 . . - yi. .. - y. } .. 4 y — 

atjk = yi 3 k. - y<3. . - y».K . - y.3k . 4 yi .. . 4 y.3 . . 4 y. . - y 


34 


Following these example deviations, one can form all the a components involved 
in the expression for uij* and thereby compute their SS's. A summary of the com- 
putational procedure is shown below. 


Space 

Source 

Spanned By 

Dimension (DF) 

SS 

Am 

— 

A 

-A 

■ » ai 

1 

SS - UKMy* . 
u — 

X A 

A main effects 

ai, .. 

1-1 

SSa JKM £ (a,V 


B main effects 

*• 

ax, . . 

~« 

a.j 

J-l 

SS* - IKM £ (a,') 2 


C main effects 

< 

ax, .. 

~c 

., a* 

K-l 

SSc IJM £ (i, c ) 2 

<. 

AB interactions 

-AB 

an , . 

"A 1 

. . t au 

(I-1)(J-1) 

SSab KM £ £ (at*/) 2 

^BC 

BC interactions 

*»c 

an, . . 

- 9C 
. . , aj« 

(J-1HK-1) 

SS«c = IM £ £ (aj* c ) 2 

^AC 

AC interactions 

* AC 

an, . 

•At 

. . , an 

(I~l )(K— 1) 

SSac - JM £ £ (a *.£)* 

^A*C 

ABC interactions 

-ABC 

axxx, . 

- A|C 

. . * a: ja 

(I-1)(J-1)(K-1) 

SSabc- M ^ £ £ (au* 1- f 


Error e 

lyiik. 

- J 

IJK(M-l) 

ss. £ £ £ £ fy:jk* - y •. jk . ) 2 

l 1 k » 

X. 

Total about 
grand mean 



1JKM-1 

SS T £ £ £ £ (y'i*» ~ y . . . .) 2 


The decomposition of Jtf into its direct sum of the above nine orthogonal subspaces 
gives rise to the following identity: 

SSt SS SSa + SSb + SSc ■* SSab + SSac + SSac + SSabc + SS, 

Once the SS's have been computed, the MS's arc readily obtained by dividing with the 
corresponding UK's. 

At this point we can answer questions concc ning the significance of each 
factor or any of their interactions. We can solve, that is, the decision theory problem. 
The various hypotheses to be tested can be set up and the tests are based on compari- 
sons of MS ratios with theoretical values of an 7 distribution. We will not go into 
further detail for this mai step since in our case we did not do any tc.'ti xg, the 
reason being we wo e only concerned with the effects in a relative sense, a fact 
which in our case could be inferred from the jcs of the SS's directly. Inste'id 
we will give a detailed description of how this method was implemented in our 
case ami present the results that we obtained . 
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The factors which we considered to affect the precision of the baselines 
were their length, the mean orbital height, the a priori weights on the station 
coordinates, the elevation cut-off angle, and the orientation of the baselines in 
the grid. Since the grid design was fixed, the same baselines would be estimated 
from each solution irrespective of the selection of the rest of the factors. Due 
to the symmetry in the grid it was possible to select a set of baselines ranging 
from 25 to 300 km in length, at 25 km increments (12 levels of variations!, for 
which we could cover all possible realizations of the rest of the factors, but mainly 
the baseline orientation variations. It can be verified from Fig. 1 that the cardinal 
directions of the grid are running almost parallel and perpendicular to the satellite 
groundtracks. We selected as our reference the side defined by stations 91 and 11, 
with respect to which we determined the orientation of the various baselines (i.e . , 
any baseline parallel to 91 -11 has bearing O', any one perpendicular to it 90', etc. i. 

To avoid overcomplicating the setup we grouped the baselines into three groups; 
baselines with bearing between 0 ; and 30', between 30' and 60', and between 60' 
and 90 : ( three levels of variations). Since we were only examining two types of 
orbits the mean orbital height had only two variation levels, 400 km and 1000 km. 

The nominal elevation cut-off angle for laser observations is usually 20'. The 
reason we were interested in an increased cut-off angle is that this would result 
in a simpler retro -reflector design with a considerable decrease in cost. Two 
levels of variations were therefore considered, the nominal 20' and the 35" option. 

The last factor to be considered, the a priori station coordinate weighting, had 
three levels of variation. A priori standard deviations of 1 m in the three 
Cartesian coordinates of all stations in the grid denote strong weighting. In the 
second variation the standard deviations are increased to 25 m denoting medium 
weighting, and the third variation is the solution with "quasi-minimum" constraints. 

The latter were applied on stations 11 (d . =0.001 m, d. = 0 . 001 m), 15 (d , = 0 . 00 1 m, 
CT ; =0.001 m>, and 95 (d , = 0. 001 m, d : = 0. 001 mi. For the above setup of five 
factors with their associated levels of variation we have a total of 432 treatments, 
half of which pertain to the low orbit, the other half to the high orbit . From the 
twelve solutions which are needed in order to cover all possible combinations of 
factor variations, we selected the standard deviations ot baselines which fell in 
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each of the treatment categories and formed the ANOVA Table 7. Only one obser- 
vation per cell (M=l) was made since we were not interested in testing any hypotheses, 
but only in ascertaining the relative importance of the five factors as reflected on 
the corresponding SS's. A summary of the setup and the notation in Table 7 is the 
following: 


Factor 


Levels 


L 

baseline length 

12 

25 km through 300 km incrementing 
by 25 km 

H 

mean orbital height 

2 

400 km (OSUL) and 1000 km (OSUH) 

W 

weighting schemes 

3 

A : a x = <7y = cr z = 1 m all stations 
B: a * = a Y = cr z = 25 m all stations 
C: "quasi-minimum" constraints (see 
text) 

E 

elevation cut-off angle 

2 

20° and 35° above the horizon 

O 

baseline orientation 

3 

0°-30° , 30°-60°, and 60°-90° (see text) 

Using the data from Table 7 we computed the SS 

's and MS’s following the standard 


computational procedure as described previously. The numerical results are 
presented in Table 8. It can be readily verified from this table that the elevation 
cut-off angle is responsible for most of the variation in the data, followed by the 
rest of the factors with considerably smaller effects. In order to determine the 
importance of the factors in the case of the two orbits (OSUL and OSTJH) independently, 
two more tests were performed. For each test only half of the data (pertaining to 
the relevant orbit) were analyzed . Numerical results are shown in Table 9 for OSUL 
and Table 10 for OSUH. These tests revealed some interesting facts for the effect 
of the orbital height on the way the rest of the factors affect the baseline precision. 

At first the dominance of the elevation cut-off angle was reconfirmed as well as the 
fact that baseline length variations are second in line as far as the precision is 
concerned. For the low orbit, however, the orientation of the baselines is more 
important than the weighting scheme, while in the case of the high orbit it is the 
one with the least effect. This is reversed for the weighting scheme which in the 
case of the high orbit is almost as important as the baseline length factor. This 
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Table 8 ANOVA Results for Test A 


ANALYSIS OF VARIANCE TEST A 

LEVELS OF FACTORS 


L 12 
H 2 
W 3 
E 2 
O 3 


GRAND MEAN 1.81250 


SOURCE OF 

SUMS OF 

DEGREES OF 

MEAN 


VARIATION 

SQUARES 

FREEDOM 

SQUARES 

£ Tu'13 si 

L 

37.42303 

11 

3.40255 

1.84 

H 

3. 4048 1 

1 

3. 484a 1 

1.87 

LH 

1. 17963 

11 

0. 10724 

0.33 

V 

S. 49847 

o 

4.24524 

2.06 

LV 

6.24431 

22 

0.28333 

0.53 

HW 

9.63083 

2 

4.81544 

2. 19 

LHW 

0.93963 

22 

0.04271 

0.21 

E 

144.90753 

1 

144.90750 

12.04 

LE 

10. 15917 

1 1 

0.92356 

0.96 

HE 

0.01815 

1 

0.01815 

0. 13 

LHE 

0.85519 

1 1 

0.07774 

0.23 

WE 

2.45097 

2 

1 . 22349 

1. 11 

LVE 

1.03069 

22 

0 . 04635 

0.22 

HWE 

4.78949 

o 

2.39475 

1.55 

LEWS 

0.34218 

22 

0.03 328 

0.20 

0 

3.63625 

o 

1.82313 

1.35 

LO 

5 . 25936 

22 

0.23908 

0.49 

HO 

2.C9083 

2 

1.04344 

1.02 

LHO 

2.6396S 

22 

0. 11999 

0.35 

VO 

2.34933 

4 

0 . 53747 

0.77 

LVO 

2.30569 

44 

0 . 05240 

0.23 

HVO 

1 . 6305 1 

4 

0.40763 

0.64 

LHVO 

1.83394 

44 

0.04168 

0.20 

EO 

3.94631 

2 

1 . 97340 

1.40 

LEO 

2.68319 

22 

0. 12219 

0.35 

HEO 

0.38332 

o 

0.44266 

0.67 

LHEO 

2.3S634 

22 

0. 10847 

0.33 

WEO 

1.16431 

4 

0.29108 

0.54 

LVEO 

1.66236 

44 

0.0C77S 

0. 19 

HWEO 

1.38329 

4 

0 . 345S2 

0.59 

LHVEO 

1 .65005 

44 

0.03750 

0. 19 

TOTAL 

269.99250 

431 




\ 
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Table 9 ANOVA Results for Test B (OSUL) 


ANALYSIS OF VARIANCE. . 

, . . .TEST B 




LEVELS OF FACTORS 





L 12 





W 3 





E 2 





0 3 





GRAND DEAN 

1.90231 




SOURCE OF 

suns of 

DEGREES OF 

MEAN 


VARIATION 

SQUARES 

FREEDOM 

SQUARES 

* EMS a 

L 

22. 18718 

11 

2.01702 

1.42 

V 

2.96037 

2 

1.4E019 

1.22 

LV 

1.56296 

22 

0.07104 

0.27 

E 

70.34116 

1 

70.84116 

3. 42 

LE 

7.58051 

11 

0.63914 

0.63 

WE 

0. 19704 

O 

0.09852 

0.31 

LWE 

0.09963 

22 

0.00453 

0.07 

0 

4.40481 

2 

2.20241 

1.43 

LO 

5.80519 

22 

0.25337 

0.51 

WO 

0.24241 

4 

0.05060 

0.20 

LVO 

0. 30093 

44 

0.0C534 

0.03 

EO 

3.68925 

2 

1 . 84463 

1 .36 

LEO 

3.78974 

22 

0. 17135 

0.41 

WEO 

0. 10296 

A 

0.02574 

0. 16 

LWEO 

0. 13370 

44 

0.00304 

0.C6 

TOTAL 

123.88884 

215 




Table 10 

ANOVA Results for Test C (OSUH ) 



ANALYSIS OF VARIANCE. 

.... TEST C 




LEVELS OF FACTORS 





L 12 

W 3 

E 2 

0 3 





GRAND MEAN 

1.72269 




SOURCE OF 

SUMS OF 

DEGREES OF 

MEAN 


VARIATION 

SQUARES 

FRLEDOII 

SQUARES 

* F31S a 

L 

16.42051 

11 

1 . 49277 

1.22 

V 

15. 16398 

2 

7.58449 

2.75 

LW 

5.62102 

22 

0.25550 

0.C1 

E 

74.08449 

1 

74.08449 

8.61 

LE 
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significant interaction of orbital height, weighting scheme and baseline orientation 
is also confirmed from the corresponding SS's (SShm and SSho ) of Table 8 for the 
complete data set analysis. The large value for SShh in this table indicates that 
the applied a priori weights on the station coordinates will produce significantly 
different results depending on the orbital height of the spacecraft. Combining this 
with the results of Tables 9 and 10 we can say that the higher the orbit the larger 
the variation in the precision of the baselines due to identical variations of a priori 
station weights. These general remarks confirmed what we intuitively expected. 

In the following we examine die data for each of the factors individually, and we 
present whenever possible the theoretical explanation for the trends exhibited in 
them. 

B. Baseline Precision Variations Due to Different A Priori Station Information 

The plots which are presented and discussed in this section are graphical 
representations of the tabular values (Table 7) of the standard deviations and 
provide a more illustrative tool for the investigation of their variations. The 
quantities denoted by SA, SB and SC in these graphs correspond to the standard 
deviations obtained from the three weighting schemes A, B and C respectively. 

The quantity SM corresponds to the arithmetic mean of the SA, SB and SC. Figs. 
5, 7 and 9 Illustrate the variation of the standard deviations for both systems, 
for the three orientation classes and for a 20° elevation cut-off angle. Figs. 6, 

8 and 10 give the same information when the elevation cut-off angle is increased 
to 35°. Figs. 11 and 12 finally were produced from the same set of data by 
averaging over the three classes of orientation. 

From Figs . 5a, 7a, 9a and 11a one can conclude that the relative influence 
of the weights remains unchanged for differently oriented baselines, at least for 
the given geometry. We can reach a similar conclusion for the case of the high 
orbit from Figs. 5b, 7b, 9b and 12a. These conclusions seem to hold true for 
either choice of the elevation cut-off angle, as can be gathered from Figs. 6a&b, 
8a&b, lOa&b, and lib and 12b. Inferences, therefore, about the influence of the 
a priori station position information can be drawn from Figs. 11 and 12 alone. 
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Concentrating on these figures we can readily verify that the standard 
deviations in most of the examined cases have a nearly linear tendency to increase 
with increased baseline length, a fact which was intuitively expected. The only 
significant exception from this rule is the case of "quasi-minimum" constraints 
for the high orbit when the elevation cutoff angle is 35°. Possible causes behind 
this result will be examined later. Relatively speaking, in the case of the low 
orbit the results of the test with "quasi -mini mum" constraints are always the best 
in terms of absolute magnitude, while for the high orbit they are of intermediate 
quality for the 20° cutoff angle test and fluctuate around 2.8 cm with no definite 
pattern when the cutoff angle is 35°. From Fig. 11a it is clear that for the low orbit 
the , 'quasi-minimum" constraints and the uniform weighting with a priori station 
standard deviations of 1 m produce almost identical results. From Fig. lib we 
can see that these two schemes produce different results, the latter being about 
10% above the first one in all cases. By changing the elevation cutoff angle, the 
number of observations was reduced equally for all three cases. A, B, and C. 

One would therefore expect an increase in the standard deviations proportional to 
the ratio /n 9 Q<> / / ngg«, with n being the total number of observations in each case. 
Such an increase would be manifested in the graphs as a shift in the "baseline 
sigma" scale and the relative location and shape of these graphs should remain 
the same. The reason they do not come out as such lies in the fact that even 
though all observations are of the same quality in terms of accuracy, they are not 
the same in terms of "geometric" quality. As it will be explained in a later 
section, the low elevation observations provide geometric strength in the solution 
which cannot possibly be compensated for by an equal number of medium or high 
elevation observations. Comparing the two figures (11a and b) one can verify that 
the loss of the low elevation observations affects the solutions with uniform weights 
more and within these the weaker estimates are obtained from the solution with 
the smaller weights. What is important from all these observations on the two 
graphs is that the internal structure of the observations plays a catalytic role on 
the performance of a given weighting scheme. The loss of observations from one 
case to the other is about 47% which warrants an increase of the standard deviations 
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by a factor of about 1.37. Even if we consider this, the resulting estimates show 
an additional 40% deterioration which is the result of the missing low elevation 
observations as explained above. 

These conslusions hold true for the high orbit also (Figs. 12a & 12b). The 
effect of the elevation cutoff angle change on the influence of the weighted constraints 
in the baseline standard deviations is even greater in this case. A 40% reduction 
in the observations implies an increase in the standard deviations by a factor of 
about 1.33; fromthe results, however, it is evident that even after we reduce the 
estimates to the same number of observations, the 35° elevation cutoff angle test 
yields poorer results by about 60% on the average. It is also interesting to note 
that strict uniform constraints produce rather optimistic standard deviations 
compared to a "quasi -minimum" constraints solution. As for the rather "irregular" 
results of the latter, in the high elevation cutoff angle test (Fig, 12b), the explana- 
tion lies in the missing observations rather than the applied constraints, the 
reason being that the constrained stations lie in the border of the grid and in the 
nominal case (20°) they collect a larger number of observations compared to the 
rest of the stations, especially at low elevations. With the cutoff angle increased 
to 35°, this advantage is lost and with it the strong "link" between the applied 
constraints and the parameters under estimation. This, in turn, produces a very 
loose system and therefore an ill-conditioned set of normal equations. An examina- 
tion of the station coordinates' quality of recovery reveals some interesting facts. 
Although the coordinates are not estimable quantities, the application of the 
constraints changes their status to what is called "conditionally estimable. " 

Since the constraints are identical for both elevation cut-off angles, a relative 
comparison of their standard deviations is not completely unjustifiable. The 
average results are shown in Table 11. 
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Table 11 


Elevation 

Cutoff 

Angle 


Average Standard Deviations (in cm) 
X Y Z 


20° 

1.0 

1.2 

1.2 

35° 

2.0 

3.5 

2.5 


It is obvious from the above table that the loss of the low elevation observations 
has affected the three coordinates differently. The standard deviations in the Y 
component increased by a factor of three compared to those for X and Z which 
only doubled. This weakness in the Y component is also evident from the compari- 
son of the correlations in the Y components among different stations in the grid 
which are approximately equal distances apart. In the 20° case these correlations 
show an increasing trend from 0. 28 for stations near the northeast side of the 
grid, to 0.88 for stations at the southwest area. In the 35° case the same trend 
is also present, but in this case the rate of change is much steeper, from about 
0.22 at the northeast side to almost 0. 99 at the southwest. The increase in the 
correlations can only be due to the loss of the low elevation observations, since 
no other parameter was varied between the two tests. As for the actual trend, 
which is present in both cases, it is probably due to the way the Y component was 
constrained: the Y coordinates of stations 11 and 15. As can be seen from Fig. 2, 
both these stations lie at the northeastemmost side of the grid and one should 
therefore expect a weaker determination of the Y components as one moves away 
from the vicinity of stations 11 and 15. The combined effect of weak determination 
and high correlations in the Y components accounts for the irregular behavior of 
the baseline standard deviations. This view is also supported from the results 
depicted in Figs. 6b, 8b, and 10b. As it is seen from the first two figures, baselines 
with orientation angles from 0° to about 60° wnh respect to the direction defined oy 
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stations 11 and 91 are much more affected than those with orientation angles 
between 60° to 90°, the reason being that due to the specific location of the grid 
in space the baselines in the second group have significantly smaller components 
along the Y-axis compared to the baselines in the first group. The errors 
therefore propagate in a more favorable way for the second group of baselines 
(Fig. 10b). 

In addition to the tests essential for the comparison of the different weighting 
schemes, some tests were conducted with different "quasi-minimum" constraints 
in connection to the rank deficiency problem m the short-arc mode. Description of 
this problem and a discussion of the results was presented in the section on rank 
deficiency and ill-conditioning in short-arc solutions. 

C. The Effect of the Elevation Cut-off Angle on the Baseline Precision 

From the ANOVA tests on both orbits (OSUL and OSUH), it is obvious that the 
elevation cut-off angle is the factor responsible for the largest variation in the baseline 
precision. Here we examine these variations in a more detailed manner, accounting 
for the fact that an increase in the cut-off elevation results in a simultaneous decrease 
in the number of observations in the problem. With the frequency of the observations 
held fixed for any choice of the cut-off elevation, the loss of observations is propor- 
tional to the percent reduction of the originally "observable" portion of a given pass. 

Based on simple geometrical relationships between the orbit and the observ- 
ing station, we can derive the following formula for the length of the observable arc, 
given the maximum elevation (E max ) that the satellite reaches with respect to the 
station and the minimum elevation (E m in) beyond which no observations are 
permissible: 

! 

g _ \ sin [(E max + E min ) + (Pi + Po)] sin [(Emax - E m in) + (Pi ~ Po)l 

sin (Emin + p o) 

where: 

S is the length of the arc in geocentric angle measure 

Po is the parallactic angle at the satellite, subtended by the geocentric radius 


53 



"1 


of the station when the satellite elevation is equal to E m i n 
Pi is the parallactic angle (as Rd) for satellite elevation E max . 

The parallactic angles are computed from the following formulae given the radius of 
the orbit a (assumed circular) and the mean earth radius a « : 

sin Po ~ ” cos E]Yifn and sin Pi — “ cos E m^y 


The arc lengths can be given in time measure also once the period of the satellite is 
computed from Kepler's third law: 


P 



GM = 398 603 x 10 9 m 3 /s 2 


The periods for die low and the high orbits in our test are for example: 
OSUL: p^! 92 min; OSUH: p 5? 105 min 


Based on the above, Table 12 gives some numerical examples for the arc length of 
overhead passes (E max = 90°) for three different choices of the minimum cut-off 
angle. In addition to the angular length and the duration of each arc, the percent 
reduction is also given for comparison purposes . 


Table 12 


Cut-off 

Angle 

Or 

bit 

OSUL 

Arc Length Duration 

OSUH 

Arc Length Duration 

0° 

20° 

35° 

40° ~ 10 min 

1 63% 

15° 75% ~ 4 min 

1 33% 

10 ■* — ~ 3 min 

60 — ~ 18 min 

i 50% 

30° 67% ~ 9 min 

i 33% 

20°- 1 ~ 6 min 


The conclusion that can be drawn on the basis of these percentages is that 
the reduction in length, hence in observations too, is very significant (33%) for the 
two cut-off angles considered in our tests . The additional implication is that we 
are not decreasing our observations uniformly (such would be the case if we had 
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decreased the frequency of the observations) but only at the beginning and the end 
of each pass . As pointed out already in the previous section, even if all our obser- 
vations are of the same precision, their geometric quality depends on the relative 
positions of the observer and the target. It is rather simple to realize that when a 
distance is to be estimated indirectly, the best measurements to do are the distances 
between its end-points and a third point on the same line. The low elevation observa- 
tions are the "third" points for our problem. It is conjectured in [Van Gelder, 1978] 
that each of these low observations contains as much geometric information as two 
observations in medium elevations. 

To get an idea of the geometric quality of these observations, we have plotted 
in Fig. 13 the average standard deviations for our baseline sample for the two 
different cut-off angles. At first glance one might think that indeed the precision is 
'■> ' by a factor of two going from the 20° (S20 curve) to the 35° (S35 curve). We 

must consider though that the loss of the observations should be accounted for first 
and is depicted by the dashed curve. So this curve gives the expected precision if in 
the 20° angle case we had decreased (uniformly) our observations to as many as we 
had in the case with 35°. On the basis of this curve, the deterioration of the results 
due to purely worse geometry is about 40%. In order to overcome this we would have 
to increase observational frequency (already 10 pps) to unrealizable rates. It is, 
therefore, recommended that we include as many low observations as possible. 

This, of course, should be further examined when systematic effects (like measure- 
ment biases, atmospheric refraction model inadequacies, etc. ) are included in 
future simulation studies. It is expected that the accuracy of these observations 
(although their precision may be the same) will in general be poor due to biases and 
unmodelled effects . One should therefore try to find the "golden cut" so that with 
proper weighting of these observations or even setting limits as to their number per 
pass, an improved geometry can be achieved with tolerable biases at the same time. 
To these direct effects of the cut-off angle we must add its interaction with 
i two other important factors. We already discussed the first one, the weighting 

( schemes. In the case where "quasi-minimum" constraints are applied, it must be 

[ established beforehand that the stations which are to be constrained have enough 
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Fig. 13 

and uniformly distributed observations on all passes in the solution. It should 
always be kept in mind that the constraints "flow" through these observations in 
order to determine the orbit, based on which the observations from the unknown 
stations determine their position. 

The other factor which seriously interacts with the cut-off angle is the 
baseline orientation relative to the satellite passes . Again the reason lies behind 
the geometric quality of the low observations . More details will be given in the 
next section where the orientation factor is examined. 
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D, Variations of the Baseline Precision Due to Different Network-Satellite 
Pass Configuration 

One of the factors which affect the precision of the recovered baselines 
in satellite ranging networks is the relative orientation of these baselines with 
respect to the sateUite pass(es). Because of the dependence of this effect on the 
adopted cut-off angle, we have already given some hints on the source of the 
problem in the previous section. These remarks, however, were based on purely 
intuitive geometric considerations. In this section the problem is examined more 
systematically, and the use of a simple example will clarify the situation and 
provide some justification for the results of the numerical tests. 

Our simple setup is shown in Fig. 14. A satellite pass lies on the plane 
defined by the axes X and Z . Disregarding time for the moment and denoting by 
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subscripts i the satellite point and j the station, the geometric distance 
between them is 

rji = [(xi - xj) 2 + (yt - yj) 2 + (z, - Zj) 2 ]^. 

Now we add one more station in the picture, station k, which also observes the 
satellite at a different epoch (close to the I th ) and measure the range r„i. From 
two sets of this type of observation we want to estimate the distance (baseline 
length) between stations J and k. For the sake of simplicity let us assume that the 
orbit is perfectly known, so the only unknowns will be the station coordinates (Xj, 
Yj, Zj) and (X*, Yk, Z k ). The design matrix of partial derivatives with respect to 
the parameters will typically look like the one below: 

Parameter: Xj Yj Zj X k Yk Z k 

r „. .In h o o o 

rjt r ji rji 

Xi-Xk Yi-Yk Z)-Z k 

^ 0 0 0 r k t ~ rit i Tki 

Examination of this matrix indicates the sensitivity of the system in each param- 
eter under estimation. We consider two different cases of baseline-pass configu- 
ration. In the first case we examine the baseline 1-2 which is perpendicular to 
the plane of the pass. Since all satellite points have zero Y coordinates, the 
derivatives 3rn/d Yi and 9 rz //9 Yz are equal to Yi/rn and Ys/t s i respectively. 
Unfortunately, the ranges do not vary too much in a short interval such as a 
ten-minute pass and therefore neither do these partials. This set-up therefore 
is very insensitive with respect to Yi and Ya , and such a solution would yield very 
poor estimates for these parameters. For this case the baseline length bia is 
simply | Yi - Ya [ , so obviously the poor results will propagate in the determination 
of bi 2 . We now examine the other extreme case where the stations 3 and 4 define 
a baseline on the plane of the pass, to this case the derivatives for Yi and Y 3 are 
zero, since Yi = 0 = Ya; so their determination is impossible. The baseline b 3 4 
is now given by 

bs4 = [(X 3 -X 4 ) b + (Z 3 - Z 4 ) 2 ]* 
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and the result therefore will be Independent of the estimates Yi and Ya . Between 
the two extreme cases discussed above, the quality of the recovered baseline 
changes from poor to best as Its orientation varies from perpendicular to parallel 
to the plane of the pass. 

One must realize that these examples are Indeed oversimplified, and In 
reality the situation is much more complicated. They are adequate enough, however, 
to illustrate how poor geometry can affect the results of equally precise observations. 
The regular grid design of the network which we investigated and the convenient 
orientation of the satellite passes with respect to its sides provided an excellent 
set-up for numerical tests. In Fig. 15 we plotted the average standard deviations 
from our standard sample of baselines, for the three different orientation classes 
which we considered as representative for this problem. We observe that in general 
the baselines which belong to the (O'-SO 0 ) and (60*- 90°) orientation classes are 
better determined than these which lie in the middle. This, of course, happens 
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because these sets of baselines are more "parallel" to the satellite passes than the 
others . A slight variation in the results for the (60* - 90*) class Is most probably 
due to a slightly unbalanced collection of observations from ascending and descending 
passes. 

Finally we must point out again that although relatively speaking the baselines 
which are nearly aligned with passes will have smaller standard deviations than die 
rest, their absolute quality will depend highly on the adopted cut-off elevation as 
already explained in the previous section. 

E. Baseline Bias Due to "Erroneous" Geopotential Model 

When we outlined the purpose of this investigation we stressed that we were 
mainly interested in the determination of the degree of dependence (or independence) 
of the baseline standard deviations on certain factors. In this sense we are unable 
toquote absolute numbers for the expected accuracy of the results except possibly 
for the component due to noise in the observations only. Average results for both 
orbits (OSUL and OSUH) are depicted in Fig. 16. Most systematic effects (e.g. , 
refraction) have been disregarded throughout the course of study. It is known, 
however, that almost all models used to correct for these systematic effects are 
imperfect and an uncertainty is always attached to their results. It is therefore 
expected that our baseline accuracy will be affected by these uncertainties and in 
fact worsened. Results on the magnitude of these components of the total variance 
were recently given by [Smith, 1978]. 

In addition to the inflation of the total baseline variance, imperfect corrections 
for systematic effects introduce biases also in the actual baseline length estimates. 

We alread* minted out that in our investigation the only case where such effects 
were of concern to us was in deriving the baseline precision from an estimation 
process based on a different geopotential model from the one used to generate 
the observations. Various cases were rerun using GEM 7 (3,1), GEM7 (8,8), and 
GEM9 (16, 16); and the baseline precisions were compared to the original tests with 
GEM7 (16,16). In all cases the results were identical down to the millimeter 
level. We noticed, however, very significant differences in the recovered 
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baseline lengths and in the height differences between stations. This was not 
further pursued since it is a major problem in itself and another study quite 
different from this one should deal with it. We do want to stress that this must 
be cleared to satisfaction before concepts such as "repeatability" are employed 
in order to determine station motions from baseline length variations. 
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7. SUMMARY: CONCLUSIONS AND RECOMMENDATIONS 


The objectives of the conducted investigation were to determine the variation 
of the baseline precision in an SRS system due to various factors such as the a priori 
weighted constraints on station positions, die increase (or decrease) of the elevation 
cut-off angle, the relative orientation of the baselines and satellite passes, and the 
accuracy of the observations. 

Through ANOVA tests for both of the orbits considered, it has been established 
that these factors produce significant variations when our accuracy goal is 1 cm. 

Aside from their direct effects, it has been shown that some of these factors interact 
with each other producing very undesirable results. The decrease of the elevation 
cut-off angle is shown to be die factor responsible for the largest variations which 
in this case are due both to the decrease in the number of observations as well as 
to the degradation of the geometric strength of the system as a whole. The effect 
of changing the a priori station information was examined in parallel with the 
problem of constraints and rank deficiency in the system. From numerical tests 
it is established that for the specific problem examined, the effective rank deficiency 
(inherent rank deficiency plus ill conditioning in the no.~rial equations) is at least 
four and on the basis of theoretical considerations at most six. The necessity 
therefore for "quasi-minimum" type of constraints is clear and so is the fact that 
their number and arrangement in the network depends highly on the geometry of 
each individual problem. As far as Bayesian (biased) estimation techniques are 
concerned, it has been explained that no meaningful results will be obtained unless 
the a priori covariance matrix for the station positions reflects reality to a high 
degree of approximation. In view of the mixed emotions in the statistical world for 
the Bayesian estimation techniques, further detailed study of the theoretical basis as 
well as the appropriateness and the consequences of such techniques applied in 
geodetic problems must be undertaken. The usefulness of estimable parametrizatton 
of our problems was pointed out at several Instances. This should be considered 
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when future software is developed or the existing software is grossly modified. 

The geometric strength of the network-sMellite passes configuration plays an 
Important role in obtaining uniform results throughout the surveyed area. In 
the design of the ground networks therefore we should always consider the kind 
of satellite coverage that the specific area has (assuming that this is dictated by 
the adopted orbit which will not change drastically in the short time frame of the 
survey). Obtaining the optimum network in this sense will not be feasible in all 
cases since other factors (such as the direction of expected motion or the actual 
location of the ares) will also impose restrictions on the design. A reasonable 
compromise should always be feasible though. In fact since most of the future 
potential users of the system have already indicated the primary areas of interest, 
it would probably be beneficial if larger simulations were conducted for all these 
areas, in these future simulation studies it is important that we include :ul known 
systematic error sources in terms of the best currently available models. For 
all study areas previous weather records must be examined and some realistic 
weather model must be developed for each one. It is rather awkward to expect 
that either all stations are visible or that none cf them is due to unfavorable 
weather conditions. The quality, quantity, and frequency of ground collected 
weather data must be established and the sensitivity of the adopted atmospheric 
refraction model must be examined in terms of the resulting biases due to residual 
refraction. A decision must be taken with respect to the type of the laser to be 
used in the operational system and Instrumental biases from laboratory calibrations 
should be Included in the future simulations. The problem of the geopotential model 
used in these simulations seems to be much more complicated than what was 
originally expected. Although the precision of the results is insensitive to any 
changes in the assumed model, the resulting biases may be orders of magnitude 
larger than the quoted standard deviations . Although this investigation was not 
concerned with the accuracy of the baseline lengths but rather with their precision, 
some of the tests which we performed using different goepotential models indicate 
that a more detailed examination of the problem is in order. The argument of 
"repeatability" which is so much used in recent publications on SRS must be 
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re-examined. We do not actually know whether slow processes (such as strain 
accumulation) will leave the local field unchanged between regular resurveys. It 
is dangerous therefore to make such an assumption because the results may show 
motions which have nothing to do with reality. In view of the use of this system by 
scientists in different disciplines, a warning must be given along with the results as 
to their suitability for the various applications. Attempts to establish strain models 
on the basis of the changes in the Cartesian coordinates of the ground reflectors 
have already iken place. In order to associate accuracy estimates with such a 
model we must first prove that strain is an estimable parametric function of the 
nonestimable coordinates. If this proves to be true, then such a treatment can be 
justified provided the coordinates are obtained from the proper estimation process 
(e.g., inner constraints least squares adjustment). 

We hope that the results obtained herein and those which are to come from 
the proposed further investigations of the system will provide sound arguments 
for the appropriateness and the capabilities of a satellite laser ranging system. 
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APPENDIX 


The following quotations serve as an introduction to the subject of this 
Appendix . 

In this name [mathematical statistics], "mathematical" seems 
to be intended to connote rational, theoretical, or perhaps mathe- 
matically advanced, to distinguish the subject from those problems 
of gathering and condensing numerical data that can be considered 
apart from the problem of inductive inference, the mathematical 
treatment of which is generally relatively trivial. The name "Sta- 
tistical inference" recognizes that the subject is concerned with 
inductive inference. [Savage, 1972] 

Subjective expectations, valuations and preferences and their 
changes from person [to person] or in the course of time can and 
should be investigated by means of "objective" statistical methods. 
Trying to use them as a basis of statistics is like trying to gauge 
a fever thermometer by means of the patient's shivers, [van Dantzig, 
1957] 

I shall call them "Bayes" probabilities because, frequency or 
not, they are the ones needed for insertion into Bayes's theorem. 
Savage argues that they are "pe rsonalistic" , that is, they are a 
property of the individual and not of society. I would dispute this 
myself, and agree with Jeffreys in saying that in scientific questions 
they are objective. They only differ between individuals because 
the individuals are differently informed; but with common knowledge 
we have common Bayesian probabilities. We can ignore this side- 
issue in the present account. [Lindley, 1958] 

I have to comment on the sentence: "It is the greatest strength 
of the Bayesian argument that it provides a formal system within 
which any inference or problem can be described". I would like to 
turn it around and say: "It is the greatest defect of the Bayesian 
argument that it provides a formal system according to which you 
can believe what you wish and, furthermore, without any data". I 
believe the search for the sort of panacea envisaged is a false one, 
which is based on a total misunderstanding of the nature of language 
and the nature of knowledge. Here again I believe some homework 
is desirable. [Kempthome, 1972] 
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Biased Linear Estimation 


The problem of estimation arises in most sciences today and although some 
still consider it as the primary and exclusive area of interest to statisticians, many 
theoretical developments can be credited to other scientists as well. The problem 
can be briefly stated as the determination under certain conditions, of quantities 
which are related through a known functional relationship to a given set of observa- 
tions. By conditions we mean a set of criteria that we establish in order to obtain 
optimum estimates in the sense implied by these criteria. A set of such criteria 
which is most often used in physical sciences and engineering is the following: 

(1) linear 

(2) unbiased, 

(3) minimum variance estimators. 

For a detailed discussion of these and other possible criteria, one may consult 
(Rao,1973). Since the choice of the criteria is more or less subjective and depends 
on the nature of the problem in hand, a good deal of controversy and confusion is 
evident from the literature whenever a comparison of different estimators is 
attempted. Most of this is due to variations of the second property— the unbiasedness. 
A number of statisticians, for instance, substitute this by "method consistency," a 
concept proposed by Fisher and Haldane. An even greater number of scientists 
follow the approach proposed by the late Professor L.J. Savage [1954 (1972 ed. ) 
and 1962], whereby they attempt to uniformly minimize the variance of the 
estimators at the expense of unbiasedness . A deeper study of the problem reveals 
the origin of the problem as being the definition of probability adopted by each of 
the parties. An elegant and extensive presentation on the four different definitions 
of probability is given in [Papoulis, 1965] . They are as follows: 

A. Axiomatic (measure theory), 

B. Relative frequency (Von Mises), 

C. Classical (favorable outcomes over total "equally likely" alternatives), 

D. Measure of belief (inductive reasoning). 

We will not argue here which of the above is most suitable as a definition, although 
we personally believe that the foundation of statistics lies in the axiomatic definition 
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rather than any empirically or intuitively conceived definitions . The purpose of 
the following sections is to present the structure of certain biased estimators 
which are often used in our areas of interest and to compare them whenever 
possible with the unbiased estimators pointing out advantages of the one over 
the other. To the best of our knowledge, the first thorough examination of 
Bayesian estimation techniques in connection with geodetic problems is [Bossier, 
1972 ]. 

The argument on which the application of biases estimation is based is that 
if we have prior information on our parameters we should use it in order to obtain 
more accurate a posteriori estimates. On the other hand, unbiasedness is essen- 
tial in our problems if we want to make correct inferences from our results. We 
should always keep in mind that geodetic problems are mainly dynamic (the earth is 
not rigid'. ). The majority of our estimates therefore are estimates of the true 
averages of the parameters over the time span of the observational data set. These 
averages change with time and the use of biased estimation techniques does not 
guarantee that the introduced biases between different solutions will be the same. 
This being the case we can readily conclude that any model for the rate of change of 
the parameters in question will be biased too. Considering that the accuracy of our 
parameters can be improved by improving the quality of our observations— which is 
possible in most cases— it seems unreasonable to seek this improvement at the 
expense of unbiased results . We can probably justify the use of biased estimation 
at preliminary stages of our research when we only want to obtain a rough picture 
of the problem with a limited number of observations. When we proceed, though, 
to explore the fine structure of the problem, such techniques should be avoided at 
all costs. 
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Best Linear Estimation 


Let us consider die model (Y, X/8, a 2 V): 

Y = X/3 + e ; € ~ N(0, a 2 V) 

and the estimable parameterlc function a T /3 to be determined . For the development 
of the theory we need not make any assumptions on the ranks of X and V. It will 
simplify the derivations though if we assume that both are of full rank. A general 
treatment of the problem is given in [Rao, 1971] . From the above set up, we may 
find a best linear unbiased estimator (BLUE) of a T /3; here, however, we are interested 
in seeing the results obtained when we drop the restriction for unbiasedness. Assume 
that the new estimator of a T /3 is b T Y. The estimation of b is based on the minimi- 
zation of the mean square error (MSE) of b T Y: 

MSE (b T Y) s E[b T Y-a T /3] 2 = a minimum (1) 

After some algebraic manipulation, we can reduce the above to the following form: 

MSE(b T Y) = a 2 fb T Vb + (b T X -a T ) (/3/o) (/3/o) T (X T b - a)] (2) 

It is obvious that we have one equation containing both unknowns /3 and O , and its 
minimization for b presupposes some knowledge for both of them or at least for 
their ratio /3/a . We have at least three choices to circumvent this problem: 

a) Use some a priori value for /3/a which we base either on prior experience (if 
any) or on what seems reasonable to us, 

b) Consider /3 as a random variable with a priori mean dispersion E[/30 T ] . Note 
that we need not know the actual distribution of /3, only the mean dispersion 
matrix is required. In this case (1) must be modified: 

MSE(b T Y| (3) = EE [(b T Y - a T /3) 2 | fi ] (Bayes' risk) (3) 

r 7 

c) Close inspection of (2) reveals the significance of the individual terms inside 

the brackets . The first term represents the variance while the second the bias 
squared. The matrix then can be interpreted as the relative weight 

that we may associate with the bias compared to the variance. One can there- 
fore select this matrix according to which of the two quantities is more important. 


A_. 
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Irrespective of which of the above or any other conceivable approach we select 
for the determination of ||| , the resulting equations are identical in 
appearance. To follow a unified approach we adopt the notation j = W 

with the proper interpretation Implied. Minimization of the MSE leads to the 
following set of equations: 

(V+XWX T )b - XWa (4) 

b - (V + XW XV XWa (5) 

b T Y = a T WX T (V + XWX 7 )' 1 Y (Bayes Linear Estimator, BLE) (6) 

b T Y = a T (X T V* X + W* )' 1 X 1 V' 1 Y (6') 

with: 

P* = (X T V* X + W' 1 )' 1 X T V 4 Y (7) 

The superscript * denotes a biased estimate. We can write die analogous expression 
for the least squares estimate (LSE) of P as: 

P = (XV 1 X)* X t Y' 1 Y (8) 

Comparing (7) with (8) we see that p*~> p as W’ 1 ■* 0. The matrix W' 1 however 
cannot be a null matrix except in special cases. If we choose for instance 
W = £rl, then W 1 = k 2 I so that as k -* 0 => W* ■* 0 in the limit, to this sense 
we can state that the LSE is the limit of the BLE . The above choice of W leads to 
a special type of BLE, introduced by Hoerl and Kennard who called them "ridge 
estimators" due to their similar mathematical structure to methods used for 
ridge analysis of second-order response surfaces [Hoerl and Kennard, 1970a and 
1970b] . These estimators will be discussed in more detail later. 

There are two points of interest that we would like to examine, namely, the 

A 

bias in P* and its mean error dispersion. Both will be compared to the LSE 
counterparts. We denote: 

T = (XVX+W V (9) 

Then (7) is written as 

P* = txVy (io) 
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To obtain the bias: 


E [$* | 3 ] « E [TXVY| 3 J 
= TXVx/3 

= T(X t v‘X + W*)3 - TW* 3 
= tt*3 - tw* 3 
= 3 - tw* 3 

In the LSE case, since 3 is unbiased, 

E [3 ] = 3 


( 11 ) 


( 12 ) 


The BLE is therefore negatively biased. It can be shown that this bias is toward 

A * 

the origin in the sense that the norm of 3* is smaller than that of 3. 

We derive now the mean error dispersion matrix for $*: 

e E t(3*- 3)(3*- 3) T ] = e !e[3*- e (3*)] [3*- e (3*)] t + 

3 Y 3 / Y Y Y 

[E ( 3 *) - 3] [e (3*) - 3l T J 

= | [O 2 (T - TW*T) + tw' 1 33 t w* T] 

= a 2 ( t - tw x T) + tw" 1 E(33 t ) w*t 

= 0 2 T - a 2 TW*T + ff*TW 4 WW*T 
= 0 2 T (13) 

where without loss of generality we have selected ff 2 W = E(33 T ). 

A 

For the LSE 3 we have: 

e ((3 - 3) (3 - 3) T 1 = a 2 (X T v* X) 1 (14) 


Recalling the definition of T we may establish the following inequality between 
(13) and (14): 

(X T V*X + W'V * (X T V*X)* (15) 

Two matrices A and B are said to fulfill A > B if A - B is non-negative definite. 
To obtain the above we must further assume that W is non-negative definite which 
is true for all three choices of W previously described. In the case that we make 

A 

a different choice of W, the mean error dispersion for 3 * is 
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( 16 ) 


D<4*) 3 E [(4*- P) (4* - /9) T ) = 0*(T - TW 4 T) + TW 4 PP r W* T 

If P = 0, then: 

D(4*) = <7 a (T-TW' l T)< ff a (XVX) 4 (17) 

For P ? 0, then, and for every choice of W, there exists a region of the parameter 
space in which (16) produces results smaller than (14). hi general, therefore, we 
can state that there is a region, including the origin (P ■ 0), for which the BLE is 
a uniformly smaller mean square error estimator of P compared to the LSE and 
another region for which die converse is true . The more general case where the 
estimate of a T /9 is not a homogeneous function of the observations Y, i. e. , it has 
the form b T Y + c with c a vector of constants, cm be found in [Rao, 1976] . 

We conclude this section giving an expression which relates the BLE with 
the LSE: 

4* = (X T V l X + W*)' 1 X T V* l Y 

= (X T V x X + wV (X T V*X) (X T V*X)* X T V* Y 
= (X T V* X + W'V (X T V* X) 4 ( 18 > 

or setting G = (X T v‘x + W 1 ) 1 X T V‘ l X, we have 

4* = g 4 (i9) 

It is obvious from (19) that the BLE can be considered as a linear transformation 
of the LSE. This is a striking similarity of this type of estimator with "shrinkage 
estimators" used to uniformly improve unbiased estimators . We point out that it 

A 

can easily be shown that the BLE "pulls" the estimate P* towards the origin as 
a whole, i.e., || 4* || s II 4 II , and not each of its components individually. 
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Ri dge Estimation 


As pointed out earlier* ridge estimation is a special type of biased estima- 
tion where V = I and W ® jr I. It is rather easy to visualise the results when 
either or none of the above is true. In [Hoerl and Kennard, 1970a] for instance, 
the case where W » K * (fl ij k 2 ] is also treated under the same name. All these 
variations can be categorized as methods using "uniform" prior distributions of 
the parameters (even though this is not explicitly stated) . We examine in what 
follows die properties of these estimators, and we compare them whenever 
possible with the BLUE. 

Under the LSE theory the expectation of die squared length (L a ) of the 
distance between the true P and its BLUE estimate P is 

E [L 2 ] = E [(P - P) T (P - P)J = a 2 tr (X T X)' 1 (20) 

We therefore obtain 

E [P T P] = P' P + a 2 tr(X T X) 1 (21) 

and assuming that the errors c are normally distributed: 

Var [L 2 ] = 20 * tr (X T Xf 2 (22) 

We are interested in the dependence of die above quantities on the condition of our 
normal equations X T X. Let the eigenvalues of X T X be X „* = Xj * X 2 * . . . 2 X p 
= X,ib > 0, where p is the number of parameters. Then: 

P P 

E[L 2 ] = O^aAi) and Var(L 2 ] = 2ff 4 £ (lAi ) 2 (23) 

On the basis of (23) we can see that die lower bounds for the average and the 
variance of L 2 are: 

o 2 /X. 1b and 2a 4 /X 2 lB (24) 

respectively. If our experiment is carefully set up bo that it fulfills the require- 
ments of a complete orthogonal design, then X T X ~ I, and we have no problem in 
obtaining stable estimates. In geodesy, however, most of our problems are non- 
linear and we very rarely have the chance to "design" the setup. These facts 
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result in extremely non-orthogonal problems which in serve ral cases produce 
"numerically singular" normals due to the high correlations among the parameters. 
The condition of a symmetric matrix is defined [Forsythe and Moler, 1967] as: 


cond(A) * || A || || A ' 1 II 

for any selected norm || • || . Specifically, for the Euclidean norm [ibid]: 
cond (A) ~ ^ 1 

The equality holds for orthogonal matrices always. The farther the condition number 
is from unity, the more ill-conditioned the matrix is. When A ItB is very small, we 
see from (24) that the distance L will tend to be large and it will vary more inten- 
sively than a slight change in our design warrants. Although our results are 
unbiased, they may be too "far" from the true values. To quote [Hoerl and Kennard, 
1970a]: "liie least squares estimate suffers from the deficiency of mathematical 
optimization techniques that give point estimates; the estimation procedure does not 
have built into it a method for portraying the sensitivity of the solution to the optimi- 
zation criterion" (minimum sum of squares of the residuals). In support of this, 
[Marquardt and Snee, 1975] go one step farther in identifying the cause of this 
insensitivity: "The 'fly in the ointment' with least squares is its requirement for 
unbiasedness." We do not comment on this since we have already pointed out the 
relevance of unbiasedness in geodetic problems. Instead we will investigate some 
interesting properties of ridge estimators . The usual form of a ridge estimate • n 
be obtained from (7) by direct substitution of V - 1 and W" 1 = k 2 1, k 2 * 0: 

0* = (X r X +k 2 I) 1 X 1 Y = KX T Y (25) 

To make matters simpler we assume thatX T X has been properly scaled so 
that it is already in correlation form and the proper transformation has been appli 
on X T Y also. Equation (25) can also be written as 

0* - [ I + k 2 (X T X)*]*/5= Z/3 (26) 

If (A t ) denote the eigen values of X T X, then the eigenvalues of R and Z, (4$) and 
( 7 ) i) respectively, are: 
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€i = l/(X,+k*) (27) 

Vt “ A,/ (A, + k*) (28) 

Since Z is symmetric positive definite, Z T Z = Z 8 , and since the eigenvalues of 
Z 2 are <T|*j ) we have: 

(0*) T (0*) = fi r Z r Zfi - 0 T Z 2 /J * T) 2 *, 0 T /3 (29) 

Now by (28) V *•» - Aj /(Ai + k a ) s 1} hence: 

II 0*1! * TU* I I’ll * II 0 II (30) 

The inequality (30) establishes the fact that ridge estimators are "shorter" in 
a global sense than their LSE counterparts. The fact that 0* depends on k 2 
mikes the residual sum of squares also a function of k 2 : 


<P*(k 2 ) = (Y -X0*) T (Y -X0*) 

« Y T Y - (0*) T X ? Y - k 2 (0*) T (0*) 


(31) 


A 

Since our criterion of optimization is to find a 0 * that gives the minimum O* 
wc see that a different minimum will be obtained for each choice of k 2 . This 
will be examined next. 
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The Ridge Trace 


Moat of the problems that we are faced with are large both in terms of die 
amoint of data pcd the number of parameters involved. In addition to this moat of 
these problems are non-linear and we very rarely see the interrelations between 
the parameters unless some geometrical approach is feasible. Ridge estimates 
claim to be the answer to such problems (at the expense of unbiasedness, of course). 

If we denote by P an estimate of P, then the residual sum of squares is: 

0 =- (Y -X/S)’ (Y-X0) 

= (Y - X/3 ) T (Y - Xfi ) + ( P - PY (X T X) (P - 9) 

- 0, u + tf(/3) (32) 

In the above 0«i B denotes the minimum obtained from the LSE theory. From the 
above we see that 0(0) is a continuous (Unction of P and the loci of O - constant 
are concentric hyperellipsoids centered at P . The continuity implies that for a 

VW 

given 0 (P ) -0o > 0 there will always be a Po which produces this 0o . Obviously 
we are interested in determining a specific Po which would fulfill «n optimality 
criterion in some sense. A natural choice is for the one that produces the minimum 
bias. We state the problem formally: 

Find a $ that minimizes P ' P under the constraint that (P - 0)’ (X T X) (P - P) 
- 0o . The solution of this problem through Lagrangian minimization yields: 

P -- (X T X + k 2 1) 1 X'Y P * (33) 

So the ridge estimator is the answer to this problem. The value of k 2 — which in the 
above is the inverse of the Lagrange multiplier for the0o constraint— is determined 
so that it is consistent with the prescribed 0o value. Actually it is simpler to 
choose a certain value for k 2 >0 and compute the 0 O value later. The previous 
approach, howc' r, sheds some light on the role of k 2 in ridge estimation. In 
(Hoerl and Kennard, 1970a] another interpretation which leads to Identical results 
is also discussed. A more general discussion which allows for rank deficient X 
and V (error variance-covariance matrix) is given in [Rao, 1971]. 
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An examination of the likelihood function for 0* reveals that the loci of 
constant likelihood are also concentric hyperellipsoids centered at 0 . Based on 
this and the similar fact that holds for the <p = constant surfaces, Hoerl and 
Kennard introduced the concept of the "ridge trace" which is the path described 
by the estimate in the likelihood space. This is practically obtained from plots 

A 

of the components of 0* for various choices of k 2 . Their justification for the use 
of this concept as a means of studying these estimators is the fact that although 

A 

"long" and "short" 0 * 's are equally likely, the 'long" ones, which will probably 
be farther away from the true value 0, may not always have equal physical meaning. 
This is where although not mathematically or statistically stated, the use of Bayes' 
approach is implied. The mathematical expression that provides useful information 
for the ridge trace is the mean square error E [L 2 (k 2 )] which is considered as the 
loss function. The following expression can be derived by use of the expectation 
operator: 

E [L 2 (k 2 )] = E [(0*- 0) T (0* -0)] 

= cr 2 h \i / (Xt + k 2 ) 2 + k 4 0 T (X T X + k 2 If 2 0 

'=i 

= Li (k 2 ) + U(k 2 ) (34) 

The first of the terms in (34) represents the total variance of the parameters, while 
the second is the square of the bias introduced when 0* is used in place of 0. To 
obtain the first result we can use the definition of Z from (26) and its associated 
eigenvalues from (28) and the fact that: 

Var [0*] = cr 2 Z(X T X) 1 Z T (35) 

In [Hoerl and Kennard, 1970a] it is shown that there exists a k 2 > 0 for 
which E [L 2 (k 2 )] < E[L 2 (0)], k 2 = 0 being the case for the LSE unbiased estimate 

A 

0. This existence theorem can be proved by examining the behavior of Li (kr ) 
and Lsfk 2 ). In the same reference the following key results are derived: 

Theorem 1: The total variance Li (k 2 ) is a continuous, monotonies lly 
decreasing function of ! : 2 . 
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Theorem 2: The squared bias La(k 2 ) is a continuous, monotonically 
increasing function of k 2 . 


t 


I 


i 

i 




On the basis of the above it is diown that: 

lim dLi(k 2 ) _ , A 2 

k^0 + d?) " " 20 (lAl > 


(36) 


lim dLafk 2 ) _ 
k 2 -»0 + d(d*) ° 


(37) 


These results are very important since they indicate that for an ill-conditioned 
system (A. lln « 0) the total variance will tend to be too large as k 2 0 + while 
the bias will be zero irrespective of the condition of X T X. As we move a little 
from the origin and k 2 > 0, we introduce a very small amount of bias (the 
derivative of L 2 is nearly flat around the origin) and at the same time we reduce 
the total variance tremendously. This is better understood from Fig. 1 which 
is reproduced from [Hoerl and Kennard, 1970a] . It is a graphical 
representation of Li , L a and their sum E[L 2 (k 2 )] . As it can be seen from this 
figure, the ridge estimate produces a uniformly smaller mean square error than 
that of the LSE, for 0 < k 2 <0.6. For k z 2 ? 0. 06, the ridge trace achieves its 
minimum. This point corresponds to the "minimum variance - minimum bias" 
estimate of /?, the one that Rao refers to as BLIMBE in [Rao, 1973], The 
determination of the k 2 value that will produce this estimate for fi is as follows. 
We treat the general case that for each component of & a different k 2 is adopted. 

Let P T AP - X T X where A is the eigenvalue matrix and P the matrix of 
the eigenvectors of X T X. Denoting ck = P0 it can be shown that starting with an 
approximation of fi and choosing each kj 2 = a 2 /a 2 , the iterative solution until 
a T Qf stabilizes will produce the desired set of k 2 's and the corresponding ridge 

a 

estimate jS*. In the case of large systems where this procedure may be very 
tedious and expensive, the direct solution using the pseudoinverse of X T as shown 
in [Rao, 1973] might be more efficient to use. Numerical examples for the 
iterative approach can be found in [Hoerl and Kennard, 1970b; Marqoardt and 
Snee, 1975]. In the second reference the results of the ridge estimator are 
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RIDOE HOMSSION 



Fig. 1 

compared to the generalized inverse estimates also. Marquardt [1970] has also 
given an excellent comparison of ridge estimators, generalized estimators and 
least squares estimators. This reference is particularly valuable, for it compares 
the ridge estimates with those obtained from an "inner constraint" adjustment; a 
procedure very popular with geodesists. The ridge estimator should by no means 
be confused with the generalized Inverse counterpart, even though the two have 
several properties in common. A numerical example solved analytically [ibid.] 
with both methods provides a number of interesting results and an illustrative 
comparison of these estimators. 
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Summary 


The concept, the foundation, and some of the most popular techniques for 
biased linear estimation were presented. The emphasis placed on ridge estimators 
is not without justification. One can very easily see their structural resemblance 
to what geodesists call "weighted constraints" adjustment. A better understanding 
of these estimators will probably help in their optimal utilization in our problems 
rather than outright rejection. It is, for Instance, true that in a short-arc solution 
even when the inherent rank deficiency is taken care of, critical geometry or too 
short passes or any combination of such factors may result in a very ill-conditioned 
and unstable system. If there is no way to obtain a linear unbiased estimate (hence 
a BLUE) or if we agree that we can tolerate a certain amount of bias, then a ridge 
estimator with controlled bias may very well be the answer to our problem . A 
last remark though should be made concerning the unbiasedness of parametric 
functions . If a minimum norm least squares g-inverse is used to obtain the BLIMBE 
of a parametric function a T /3, then the resulting estimate will be unbiased if we 
erroneously assumed that it did not admit a LUE initially and if a T /J is estimable. 

For the ridge estimator this is not true. In this case, therefore, all parametric 
function will be assigned biased estimates irrespective of their status under the 
classical least squares theory. 

The epilog: 

"The science of statistics is essentially a branch of Applied 
Mathematics, and may be regarded as mathematics applied 
to observational data." [Fisher, 1925] 
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